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We study Andreev billiards of box and disk geometries by matching the wave functions at the 
interface of the normal and the superconducting region using the exact solutions of the Bogoliubov- 
de Gennes equation. The mismatch in the Fermi wavenumbers and the effective masses of the 
normal system and the superconductor, as well as the tunnel barrier at the interface are taken 
into account. A Weyl formula (for the smooth part of the counting function of the energy levels) is 
derived. The exact quantum mechanical calculations show equally spaced singularities in the density 
of states. Based on the Bohr-Sommerfeld quantization rule a semiclassical theory is proposed to 
understand these singularities. For disk geometries two kinds of states can be distinguished: states 
either contribute through whispering gallery modes or are Andreev states strongly coupled to the 
superconductor. Controlled by two relevant material parameters, three kinds of energy spectra exist 
in disk geometry. The first is dominated by Andreev reflections, the second, by normal reflections 
in an annular disk geometry. In the third case the coherence length is much larger than the radius 
of the superconducting region, and the spectrum is identical to that of a full disk geometry. 
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I. INTRODUCTION 

Recent technological advances in manufacturing almost ballistic semiconductors of mesoscopic size (two dimensional 
electron gas, 2DEG) coupled to a superconductor has initiated a growing interest in considering the phase-coherent 
transport and the excitation spectrum in hybrid superconducting nanostructures. In these systems the electron can 
coherently evolve into a hole (and vice versa) at [the interface between the semiconductor and the superconductor. 
This mechanism had been discovered by Andreevlil. Classically, the particle momentum is conserved and the hole is 
reflected back in the same direction as the incoming electron. This scattering process is known as |fctro-rcflcction or 
Andreev reflection. The overview of the recent progress in this field can be found in several worksnTJ. The Andreev 
scattering resulting in a discrete spectrum of single-particle,-excitations of a layer of normal metal in contact with 
superconductors on both sides was first discussed by Andreevtl Based on the Bogoliubov-de Gennes equationcl (BdG) 
the excitation spectrum (Andreev states^ of a normal metal (N) attached to a superconductor (S) was first considered 
by P. G. de Gennes and D. Saint-JamesQ. A ballistic normal metal weakly coupled to a superpppductor is commonly 
called an Andreev billiard. Such systems have been extensively studied over the past ten yearsoTJ. The bound states 
were studied e.g. for SNS junctionaiTlia. A semi=infinite N region in contact with a semi-infinite S region in a strong 
magnetic field was investigated by Hoppe et al.Ej 

The excitation spectrum of Andreev billiards depends on whether the N region is classically chaotic or 
integrableauJlla. It has been shown that a gap opens in the spectrum for chaotic billiards, while for integrable 
billiards the spectrum is most likely to be gaplessu. In these works the spectrum was calculated only close to the 
Fermi level. However, the spectrum shows interesting features throughout the entire energy range below the super- 
conducting gap. The excitation spectrum is determined in this paper for two specific Andreev billiards, namely the 
NS box and disk systems. In NS box systems a rectangular N region is in contact with the superconductor, while 
the NS disk system consists of a circular S region surrounded concentrically by a circular N region. Both systems are 
integrable and, in accordance with earlier findings, the density of states (DOS) is gapless. We shall show, however, 
that for higher energies the DOS has singularities located at equal distances from each other. 

One of the central issues in this paper is the investigation of these singularities in the DOS. We calculate the spectra 
of these systems exactly by solving the BdG equations, taking into account the non-perfect interface (mismatch in 
the Fermi wave numbers and the effective masses of the Majspal metal and the superconductor, and tunnel barrier 
at the interface). For a narrow NS junction it is justifiedH'B'El that the pairing potential can be approximated by a 
step- like function (zero in the N region and some constant value in the S region). Note that the two systems have the 
common feature that the BdG equation is separable and one of the separated wave function sets ('transverse modes') 
is the same in the N and the S region. The exact quantization condition can be very simply expressed in terms of a 
phase $ m (e) which is shown to be related to the classical action of an electron moving inside the N region between 
two successive bounces at the NS interface (see Eq. (prf)). 
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Several authorsS ! li3lij EJ have already derived the Bohr-Sommerfeld approximation for the density of states. The 
DOS can be written in terms of the probability distribution P(s) of the classical path length s between two subsequent 
bounces of the electron at the NS interface. Starting from our exact quantizations-condition expressed in terms of 
the phase 3> m (e) we re-derive the commonly used Bohr-Sommerfeld approximatiorootj of the DOS in the case of 
NS box and disk systems. Moreover, we give an analytical expression for the probability P{s) in terms of the phase 
$m(e). 

It is shown that, in the framework of Bohr-Sommerfeld quantization, the singularity of the DOS is a direct conse- 
quence of the singular behavior of P(s). Using the expression for P(s) we can reproduce the counting function very 
accurately in the entire range of the spectrum and locate the singularities of the DOS for NS box and disk systems. 
The probability P(s) depends on the geometry of the N region and the location of the NS interface. We shall show that 
from a semiclassical point of view some special trajectories of the electron play crucial roles in the singular character 
of the DOS. Based on our semiclassical analysis a simple formula is given for the location of these singularities. This 
formula gives an excellent agreement with our numerically exact calculation of the spectrum. Furthermore, it predicts 
very accurately the number of the singularities and their location in the system studied by de Gennes et al.Q We. 
propose an explanation, based on our theory about the singularities, for the pronounced peaks found by Ihra et al.tB 
(for NS systems) in their numerical calculationSj-of the DOS. We believe that the singularities found by Lodder and 
NazarovE3 (for Andreev billiards), and Sipr et alEd (for SNS systems) are related to some special classical trajectories 
of the electron for which the probability P(s) is singular. 

The other issue studied in this paper is the Weyl formula for our NS systems. In the quantization of a normal 
billiard, the counting function N(E) gives the number of levels wJioap energy is less than or equal to E. The smooth 
part of the counting function N(E) is given by the Weyl formulatjH, which has already been calculated for billiards 
of various shapes. We present a new Weyl formula for our NS box and disk systems. As expected, the counting 
function obtained from the exact (numerical) calculations 'oscillates' around the curve given by our Weyl formula. 

Our exact quantum mechanical results reveal a more complex energy level structure for NS disk systems than for box 
systems. Semiclassically, two kinds of modes exist in the disk geometry. First, the electron can hit the superconductor 
undergoing an Andreev reflection (hereafter called Andreev states), and second, the electron trajectory may not reach 
the superconductor (hereafter called whispering gallery modes) . The two modes play crucial roles in the energy levels 
of the system. Contrary to NS box systems where the DOS is proportional to the energy close the Fermi level, in 
the case of NS disk systems it is constant due to the whispering gallerv_modes. In studying the spectrum of the disk 
geometry, we have been motivated by Bruder and Imry's recent worktj. Based on whispering gallery modes, they 
proposed a physical picture to interpret the significant paramagnetic effect observed in recent experimental! Thus, 
a careful analysis of Andereev and whispering gallery states as presented in this-paper and the inclusion of the flux 
induced phase may shed further light on the issue raised by Bruder and ImrycJ. Our work on this problem is in 
progress. 

We also found that, depending on the parameters of the NS disk system, the spectra can belong to cither of three 
classes: (a) 'mixed phase', in which Andreev states coexist with whispering gallery modes and the energy dependent 
coherence length £(e) in the superconductor is much less than the radius Rs of the S region, while the DOS is singular; 
(b) the opposite case, i.e. £(e) 3> Rs when the spectrum is identical to that of a normal billiard whose radius is that of 
the outer circle; (c) when fcp£(0) is order of one (here k-p is the Fermi wave number) or the NS interface is not perfect, 
and so the energy spectrum corresponds to that of an annular billiard (the circular S region is cut out). It turns out 
that there are two relevant parameters to make this classification, fcp^(e) and k-pR§. Using these parameters we sketch 
a so-called 'phase diagram' for the three above classes. It was reported in the works on the paramagnetic effected 
that the coherence length exceeds the size of the superconducting region. Thus, according to our phase diagram it 
is possible that the whole NS disk system behaves as a full normal disk. Then a larger paramagnetic effect can be 
expected than that predicted by Bruder and Imry who included only the whispering gallery states. However, further 
work is necessary to clarify this scenario. 

The text is organized as follows. In Sec. |0] a quantization condition (secular equation) is derived from the matching 
conditions of the wave functions at the interface of the NS systems. Owing to the symmetry of the BdG equation, 
the secular equation can be expressed in a very compact form valid both for NS box and disk systems. In Sec. Ill 



the Weyl formula is derived for NS box and disk systems and compared with the exact results. Our theory for the 
singularities in the DOS is presented in Sec. Ew We sketch the so-called 'phase diagram' for the NS disk systems in 
Sec. |v|, and the conclusions are given in Sec. [V| . 
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II. SECULAR EQUATION FOR BOX AND DISK 



In this section we derive a secular equation determining the energy levels for box and disk geometries shown in Fig. 
|l|. It is possible to treat both problems in a common framework by introducing the generalized coordinates 



(xi,x 2 ) 



(x, y), for box, 
(r, ip), for disk. 



(1) 



The normal region is in contact with a superconducting region at x\ = xns (where xns = for box and xjsis = Rs 
for disk). 





FIG. 1. The normal system (N) is in contact with the superconducting region (S). The two geometrical arrangements are 
(a) box geometry of side lengths d and W' and (b) disk geometry (two concentric circles of radii Rs and i?jv). For the box, the 
generalized coordinates are (xi,X2) = (x,y), while for the disk (x\,X2) = (r,cp) centered at the origin of the superconducting 
circle. 



At the interface x\ = xns the tunnel barrier is modeled by a Dirac delta potential V(xi,x 2 ) — UqS(xi — xns); 
as in Ref. ^5|. The superconducting pairing potential is a constant Ao in the S region and zero in the N region. 
The self-consistency of the pairing potential is not taken into account just like in Ref. ^5|. For a planar geometry 
the difference between the self-consistent gap and the step function is quite small as it was shown by McMillan and 
KieselmannEa. Similar results has been found by Plehn et al. in Ref. |27] for superconducting multilayers. For disk 
geometry the step function approximation is still reliable provided Rs is larger then the coherence length. HoweAfer, 
the Fermi energies (i.e. the energy differences between the band edges of the N/S region and the chemical potentialES) 

and the effective masses in the N and S regions are assumed to be different, i.e. E^ ^ E^ and ton Wg. 
The NS system is described by the BdG equation: 



H A 
A* -H 



(2) 



where VP is a two-component wave function, Ho — p /2m — /i, and the Fermi energy is fi = E^ N) in the N region and 

/i = Ep^ in the S regioa. The energy levels of the Andreev states are the positive eigenvalues e of the Bogoliubov-de 
Gennes (BdG) equationtl. In what follows, we consider the energy spectrum below the superconducting gap, e < Ao. 
The two-component wave functions 5* in the N and S regions can be chosen in the form 



am\m' h \xi) 



Xm(x 2 ), 



(3) 



*L S) (^2) 



{xi 



^(*l) 



Xm(x 2 ), 



(4) 
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where the 'transverse' wave function of the rnth mode is 



j sm(mny/W), for box, 

Xm{x 2 ) = < . (5) 
( e im¥ \ for disk. 

The 'transverse' wave functions in the N and S regions are the same, and it is possible to separate the variables X\, x 2 
in the BdG equation, i.e. Xmi^'i) depends only on the coordinate x 2 . 

The wave function tfm' e ^ {x\) is the electron-like component of the eigenspinor of the BdG equation in the N region. 
It satisfies the one-dimensional Schrodinger equation obtained by separating Xm{x2) in the BdG equation: 

J sin km \d- x), for box, 
^ {Xl) ~ \ Jm(k e r) - ^^ Y m (k e r), for disk, 

where the energy dependence of the wave number km\s) for box and k e (e) for disk are given by 



jfc<>>(e) = 4 N) \A + £ /4 N) - {™/kf ] W), for box, (7) 



k e (e) = 4 N) V 1 + £/#F N) , for disk, (8) 

kp and are the Fermi wave number and the Fermi energy in the N region, respectively, while J m {x) and Y m (x) 
are the Bessel and Neumann functions of order m. As we shall see below it is convenient to use real wave functions for 
<fini' e \xi). The wave functions have been chosen such that they satisfy the Dirichlet boundary conditions at x = d 
for box and r = for disk. Due to the symmetry of the BdG equation the hole-like component of the eigenspinor 
of the BdG equation in the N region is 

V £> h Hs,x 1 )=^(-e,x 1 ), (9) 

where the explicit dependence of the energy e in the wave functions is emphasized for clarity. 
The wave function in the S region is 

^ s , e)(£ ^ ) = jexpH^*), for box, 



J m (q e r), for disk, 



where 



2L e) (e) = 4 S) V 1 + V - (W4 S W, box, (11) 



q e {e) = 4 b Vl + 7?, for disk, (12) 
k-p and Ep are the Fermi wave number and the Fermi energy in the S region, respectively, and 



n=yje*-&%/E!$>. (13) 

Note that for box ifm' e \x) — ► as x — * — oo satisfying the boundary condition at — oo. On the other hand, only the 
Bessel function can be chosen for disks, since the Neumann function is singular at the origin. Again, owing to the 
symmetry of the BdG equation, the hole-like component of the BdG spinor in the S region is 

^ h \e,x 1 ) = [^{-E,x 1 )]\ (14) 

Finally, in Eq. (g) 

7e = A /( e -^ 2 -A§), (15) 

lh = lt- (16) 
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The coefficients dm , a m ! \ Cm\ c^J' in Eqs. are determined from the boundary conditions at the interface of 

the NS system (see, e.g., Ref. |28|): 



tf(N)| 



d 



$(N) _ HM$(S) 
m s m 



*(S)| 

2?7lN 



(17) 



The matching conditions yield the following secular equation for the eigenvalues e of the NS system for fixed mode 
index m: 



(JV,e) 










(N,h) 
(fin 




(jV,e)~ 
(fin 


1 




7e 





(N,h) 
fin 


1 



(S,e) 

(S.e) 

(fm 



(S.h) 

lh(Pm 
(S,h) 
(Pm 



7e Zip m ' + 



Zp^> + 



P 



(S,< 
¥>rra 

(S,e) 



7h -Z^™ ' + 



0, 



(18) 



where Z = (2j7iN/?i 2 ) ?7o is the normalized barrier strength, and the prime stands for the derivative with respect to 
x\. All the functions are evaluated at x\ = xns- 

Using the fact that the wave functions pm ,e ^ given in Eq. (Q) are real functions and the symmetry relations between 
the electron-like and hole- like component of the BdG eigenspinor given by Eqs. (^|) and (fTi]), the above determinant 
can be simplified. One can show that the secular equation reduces to 



Im 



{ 7c D^( £ )D( n h) (e)}=0, 



where 



(JV,e) 
<Pm 

(Pm 



Zip. 



(S,e) 
Pm 

(S,e) m N 

ms 



(19) 



(20) 



is a 2x2 determinant, and Dm\e) = Dm\—e) ■ The energy levels of the NS systems can be found by solving 

the secular equation ( |l9| ) for s at a given quantum number m. The secular equation ( |l9] ) is exact in the sense that 
the usual Andreev approximation is not assumedo. The Andreev approximation is valid only for-Ao/^F <C 1 and 
quasi-particles whose incident / reflected directions are approximately perpendicular to the interfaces 



A. Secular equation for box 



To find the energy levels for a box we shall give an explicit form of the secular equation (19). Inserting the wave 
functions given in Eq. (|J) into Eq. (|20| ) we obtain 



Hence, the secular equation given in Eq. ( |l9| ) becomes 



Im \ 7e ( Z - 



ms 



klV 



■i^qW+k^cotk^d" 
ms j 



= 0, 



(21) 



(22) 



where km\e) = k)n {— e) and qm'(s) = Qm(— e) are the wave numbers of the hole-like quasi-particles in the N 
and S regions, respectively. For a given quantum number m the energy levels can be found by solving the secular 
equation. Note that the case sinfcm^ = corresponds to the secular equation of the entire box with Dirichlet 
boundary conditions. Thus, it cannot be zero for the NS system. Semiclassically, small wave numbers k!^ and 
in the N region correspond-ito quasi-particles incident at grazing angles on the NS interface. In this case the Andreev 
approximation is not validcl. However, the secular equation (22) is still exact for the energy levels for box geometry. 



(h) 
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B. Secular equation for disk 



We now derive the explicit form of the secular equation for disk geometry. Inserting the wave functions given in 
Eq. (|) into Eq. $2(\) we obtain 



J m {q e Rs) 



felW F m(fcei?s) ZJ m {q e R s ) 



(23) 



where the primes denote the derivatives of the Bessel functions with respect to their arguments. For a given angular 
momentum quantum number m the energy levels are the solutions of the secular equation given in Eq. (|l9|) . 



III. WEYL FORMULA FOR NS SYSTEMS 



For normal systems the counting function N(E) is the number of states whose energy is less than or equal to E. 
The derivative of the counting function with respect to the energy is the density of states. The smooth part of the 
counting function N{E) for a cavity was first derived by WeylEll (for more details see Ref. ^||). The leading term 
of N{E) is given by the integral of Q(E — H c \(p, r)) over the phase space divided by h 2 , where 8 is the Heaviside 
function and H c \ is the Hamiltonian of the corresponding classical system. In two dimensions, and excluding the 
factor 2 for the spin, the smooth part of counting function N(E) is given by 



H cl (p,r)<E 



d 2 pd 2 rQ(E- H cl (p,r)). 



(24) 



For Andreev states of energy e less than Ao it is not possible to define the classical Hamiltonian H c \, therefore the 
smooth part of the counting function N(e) cannot.he derived from Eq. (p4]). As an alternative method for calculating 
the DOS one may start from the secular equationoE]. In this section we derive a Wcyl formula for the NS systems 
shown in Fig. [l] using the secular equation ([il)|). 

For a given quantum number m let us introduce the eigenphase <fr m (e) for the NS system: 



(25) 



Here we have assumed that Dm\e) ^ 0, and in the following it will also be supposed that Dm\e) ^ 0. In Sec. |v| 



the case Dm \e)Dm l {s) = will be discussed. It is also assumed that the eigenphases $ m (e) are all monotonic and 
increasing with e. The secular equation (|l^) can be further simplified: 

$ m (e) - $ m (-e) - 2 arccos = 2mr, (26) 



i0> 



where n = 0, ±1, ±2, • • • . Equation (|26|) is a very simple quantization condition for the NS system. The solutions of 
this equation give the energy spectrum e mn below the gap. The above-given quantization condition is a convenient 
starting point to calculate the DOS and the smooth part of the counting function. However, for numerical purposes 
Eq. ( |l9| ) may be more suitable. Using ( |2^ ) the density of states p(e) = J2m=i J2^=-oo — £ mn) for the NS system 
becomes 



Pie) 



M(e) 

E 



dF m {e) 



de 



J2 5{F m {e)-n), 



(27) 



where 2irF m (e) — $> m (e) — $ m (— e) — 2 arccos(e/A ), and M(e) is the number of 'transverse modes' depending on e. 
Since the $ m (e) are all monotonic and increasing with e, the derivati ves dF m (e) / de are positive, and so the modulus 
sign is superfluous in d27n . Applying the Poisson summation formular 1 lr 2 l to the second sum one finds 



M(e) 



dF m (e) 
de 



2^ cos(2nkF m (e)) 



fc=i 



(28) 



G 



The DOS can be separated into two parts: the smooth part, i.e. the k = term, and the oscillating part, which comes 
from the terms k > 1 in Eq. (|28|). Here we consider only the smooth part. Then the Weyl formula, i.e. the smooth 
part of the counting function for NS system can be obtained from N(e) = p(e') de', and one finds 

1 M(e) /ll \ 

= 2; £ [*m(e)-* m (-e)]+M(e) ( - -- arccos^-j . (29) 

m— 1 ^ ' 

This is our Weyl formula for the NS systems. A similar treatment has been mentioned by Schomerus and BeenakkerEl 
in deriving the Bohr-Sommerfeld approximation of the DOS for Andreev's billiards. The number of transverse 
modes M(e) is a discontinuous function of e, therefore in Eq. (|2^) the sum has to be replaced by an integral to get 
the smoothed version of the counting functional One can see that the Weyl formula is expressed in terms of the 
eigenphases $ m (e) defined in Eq. (po|). Note that our Weyl formula is only valid for those systems in which a common 
set of the 'transverse wave' functions can be separated from the wave functions of the N and S regions. Lifting this 
condition is a possible extension of this problem. 



A. Weyl formula for box (no mismatch) 

We now consider the Weyl formula for a perfect interface, i.e. for ru = 1, T v = 1, Z = 0. For simplicity, in the case 
of a perfect interface we shall use the notation fcp = kL = kL and E-p = E^ — E F S \ From Eqs. (|2^) and J25|), 
and in Andreev approximation (Aq/E f <C 1, i.e. km ~ Qm except when they appear in an exponent), one finds 

* TO (e) = 2k$d. (30) 
Thus, from Eq. ( p6[ ) the quantization condition becomes 

(k$ - d = nir + arccos(e/A ). (31) 

Note that this result can also be derived from the Bohr-Sommerfeld quantization rule. In fact, h$ m (s) is the classical 
action for the electron moving along the x direction between the superconductor and the wall at x = d. 

The exact counting function N(e) obtained (numerically) from Eq. ( |l9| ) and the Weyl formula given in Eq. ( p9| ) are 
plotted in Fig. |^. For parameters see the figure caption. 

800 | . 1 . 1 . 1 . 1 . 1 



Exact N{£) 

Weyl 



1 1 1 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 

e/A„ 

FIG. 2. The counting function N(e) obtained from the exact quantum mechanical calculations (solid line) and the Weyl 
formula given in ( p9| ) (dashed line) as functions of e/Ao for the box geometry. The parameters are d/W = 3, kpW/n = 87.9 
and Ao/Ef = 0.02. 

One can see that the exact counting function oscillates around the curve obtained from our Weyl formula. 

To find an analytical expression for the Weyl formula, the sum in Eq. ( p9|) is replaced by an integral and we get 

N(e) = - p {N) E F g(e/E F ) + M h [1/2 - 1/tt arccos(e/A )] , (32) 

7T 
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where 



g(x) = y / 2x(l - x) + (l+x) arcsin - x)/(l + x) - vr/2 (1 - x), 



(33) 



and p( N ) = 2m /h 2 A /{An) is the density of states for the isolated N region of area A = Wd. The number of open 
channels for the hole-like quasiparticles is Mh = [k^W/n yl — e/E-p], where [• • •] stands for the integer part. For 
i<Cl (typically x < 0.1) g(x) « ttx. Thus, the leading term of the smooth part of the counting function of Andreev 
states is N(e) ~ 2e p^ N \ Electron-like and hole-like quasiparticles make equal contributions to the Weyl formula $33), 
thereby the factor 2. The exact counting function N(e) and the analytical form of the Weyl formula given in Eq. ( |32| ) 
are plotted in Fig. ^. It is seen from the figure that the analytical expression (dashed curve) deviates only for energies 
close to the gap. 




FIG. 3. The exact counting function N(e) (solid line) and the analytical expression of the Weyl formula given in 
(dashed line) as functions of e/Aq for the box geometry. The parameters are the same as in Fig. 0. 



B. Weyl formula for disk (no mismatch) 



In this section an analytical expression of the W eyl fo rmula is derived for disk geometry. However, most of the 



results obtained in this section will be used in Sec. IV B 



too. 

1, 



Again, we shall take the case of perfect NS interface, 
namely no mismatch and tunnel barrier are assumed (r* = 1, r v = 1, Z = 0). To find an analytical expression for 

the eigenphase $ m (e) defined in Eq. ( p^ ) we have to approximate the determinant Dm (e) in (|2^) . The details of the 
approximations are presented in Appendix [a|. 

To summarize, according to the approximations of the determinants appearing in the secular equation (fl9|), three 
ranges of to > can be distinguished: for Andreev states to < k^Rs — \/khRs, for whispering gallery states m > 
k e Rs + ^keRs, and for intermediate states k^Rs — \/khRs < to < k e R$ + \fk e R$. Owing to the degeneracy of the 
±771 states the three ranges for 777 > and to < are symmetrically located with respect to the state m = 0. This 
structure of the energy levels will be called a 'mixed phase' (MP). The exact energy levels, e nm obtained (numerically) 
from df) are shown as functions of the angular quantum number to in Fig. ||. For each to the solutions are labeled by 
the quantum number 77 (radial quantum number). Since the pairs ±m are degenerate, only m > states are plotted 
in the figure. 
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100 200 300 400 500 

m 

FIG. 4. The exact energy levels, e nm (dots) in units of Ao as a function of the angular momentum quantum number m. 
There is no mismatch, and the potential of the tunnel barrier is zero. The parameters are Rs/Rn = 0.7, Ao/Ef = 0.1, and 
fcpi?s = 350. The intermediate m states are located between the two dashed lines. The Andreev states are to the left of the 
left dashed line, and the whispering gallery modes are to the right of the right one. 



The two dashed lines in Fig. |^ separate the three types of states characterized by m. One can see that the 
intermediate states indeed occupy a narrow range in m compared to Andreev and whispering gallery states. It is also 
clear from the figure that the energy levels for Andreev states 'oscillate' with increasing amp litude as m increases 
towards the onset of the intermediate states (left dashed line in the figure) . Inserting Eq. ( |A10| ) into the quantization 
condition given by Eq. (26) one can approximately calculate the energy levels for Andreev states. Plotting these 

-levels 'oscillate' around these 
]. For whispering gallery states 



energy levels in Fig. [4| (not shown in the figure) we found that the exact energi 
approximate ones. Similar behavior of the energy levels has been found by Sipr et al.E 
the approximate secular equation is given by Eq. (All), and its solutions coincide very accurately with exact energy 

^.obtained 
3. Indeed, 



levels. These approximate energy levels for Andreev and whispering gallery states are those which can 
from the semiclassical quantization of the Schrodinger equ atio n describing the radial motion of the electro 
one can check that the first two terms for <S> m (e) in Eq. (A9) multiplied by h give the radial action for the electron 
moving between the two circles. The last three terms are constant, i.e. independent of e therefore they do not play any 
role in the dynamics of electron. These trajectories of the electrons correspond to Andreev states. The semiclassical 
orbits for Andreev and whispering gallery states are shown in Fig. [| 





FIG. 5. The corresponding semiclassical orbits for (a) Andreev states and (b) whispering gallery states. 



Inserting ( A10) into the Weyl formula ( y,9\j one can determine the contributions of the Andreev states to the smooth 
part of the counting function: 

N A s{e) = ^ l#m( k eRK)-tim(keRs)-ti m (khRN)+#m(k h Rs)]+2M A s(~-^ arccos^-) , (34) 



|m|<M AS 



where Mas = [khRs — \/khRs ] is the highest m for the Andreev states. 

For whispering gallery modes the eigenphase ^ m {e) in fl25| ) cannot be determined because Dm (e) is approximately 
zero. In t his c ase the contribution to the smooth part of the counting function will be calculated from the secular 
equation (All) for whispering gallery states. It is well known that the zeros of the Bessel functions can be determined 
to a good app roximation by applying Debye's asymptotic expansion given in Eq. (Al). The solution of the secular 



equation (All) for whispering gallery states is then equivalent to the solution of the following pair of equations for e: 







TT 4 



(35) 
(36) 



where ni,ri2 = 0, ±1, ±2, • • •, and , & m (x) is given by Eq. (A5). From these equations one can find the density of states 
of whispering gallery states much in the same way as from Eq. (p^) for the Andreev states: 



m=M„ 



ds 



k h R N 

E 

ra=M w . 



rf(e) 



ds 



(37) 



where M wgs = [fc e i?s + \/k e Rs] is the smallest angular momentum quantum number m for whispering gallery states, 
and the factor 2 corresponds to the ±m degeneracy. The contribution of whispering gallery states to the smooth part 
of the counting function can be obtained from N wgs (e) = J Q Pwgs(e') de' which yields 



2 fee-RN 2 khR N 

Nw gs (e) = - V i? m (fc e i? N ) - - V & m (k h R v ). 

TT ■ ' 7T ' — ' 



(38) 



=M W 



=M W 



The negative sign in front of the second term comes from the fact that the derivative of # TO (fc/i-RN) with respect to s 
is negative. 

The counting function N(s) for the NS disk system is the sum of the contribution of Andreev states (Nas( £ )) an d 
that of whispering gallery states (N wgs (s)). We now neglect the contribution of the intermediate states. To find 
an analytical expression for N(s), the summation over m is replaced by an integral in Eqs. (|34| ) and (|38| ) and for 



k F R s > 1 Mas « M w 



khRs is taken. The last approximation corresponds to replacing the intermediate states 



by whispering gallery states. This is quite a good approximation, since the number of intermediate states is much 
smaller than that of the Andreev and whispering gallery states. The resulting integrals can be performed analytically. 
After some algebra we obtain 



N(e) 



k w R 



' S 5 (^) + 2 fcF i ?s Jl--l---arccos- 



1 1 



where 



3 i j i _ x 

g(x) = —\/2x(l — x) H — (1 + x) arcsin \ (1 — x) 

7T 7T V 1 + X 




1 



(39) 



(40) 



This is our Weyl formula for NS concentric disk systems. Usually, e/E-p <C 1 and for i« 1, g(x) — x in leading order, 
therefore the first two terms in Eq. (|39|) yield 2p^s, where = 2m/h 2 A/(4ir) is the density of states for a normal 



annular billiard of area A = (i$ - Rj) tt. The factor 2 comes from the electron/hole contributions for NS systems. 



The higher order corrections in g{x) become more relevant with increasing s. The last term in (|39j) is related to the 
phase shift due to the Andreev reflection at the NS interface. It is important to stress that the Weyl formula given in 
Eq. (^) was derived for the case Rs ^> £(e), i.e. when Andereev states and whispering gallery states coexist (mixed 
phase). In the opposite case, as we shall see in Sec. 0, no Andreev states exist. 

In Fig. H the exact counting function and the Weyl formula given by Eq. ( |39"| ) are plotted as functions of s/Aq. For 
parameters see the figure caption. 
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FIG. 6. The exact counting function N(s) obtained from Eqs. (^) (solid line) and N(e) given by (^) (dashed line) as 
functions of e/Ao for the disk geometry. There is no mismatch and the potential of the tunnel barrier is zero. The parameters 
are: Rs/Rn = 2/7, A /E F = 0.05 and k F Rs = 100. 

The singularities of the DOS (derivative of the counting function with respect to e) will be discussed in the next 
section. 



IV. SINGULARITIES IN DOS 

The exact counting function, N(e) has cusps at some energies as shown Figs. and ^| for box and disk geometries, 
and perfect interface. This implies that the DOS of the Andreev states is discontinuous here. Moreover, these cusps, 
are at equal distances. Similar singularities of the DOS for NS systems have already been found by de Gennes et al.0 
as well aSpin the numerical works of Richter et al.c3 (for NS systems), Nazarov et al.ta (for Andreev billiard), and 
Sipr et al.EJ (for SNS systems). In this section we investigate the singularities in the DOS for the NS systems shown 
in Fig. |. 

Several authorsalljTJ have already derived the Bohr-Sommerfeld approximation to the density of states, 



poo °° / 

Pm{e) = M dsP{s)Y^8[e- 



1 ^ nhvp 
2 



(41) 



where M is the number of modes in the superconducting leads connected to the normal billiard, and P(s) is the 
classical probability that an electron entering the billiard exits after a path length s. Starting from our quantization 
rule given in Eq. ( p6| ) we re-derive Eq. (El]) in this section, and give an explicit expression of the probability P(s) 
in terms of $ m (e). We shall show that P(s) is singular at some s resulting in the singularities of the DOS. These 
singularities correspond to the cusps in the integrated DOS N(e). It turns out that the simple Bohr-Sommerfeld 
approximation to the DOS is a good approach to understand the singularities in the DOS provided that P(s) is 
approximated correctly for small s. 

The solution of Eq. (|2q) gives the discrete energy levels for the NS systems. Since e <C E-p, one can Taylor expand 
the LHS of Eq. (|26|) in terms of e and find 



(42) 



where "^(O) denotes the derivative of $ m (£) with respect to e, evaluated at the Fermi energy, i.e. e = 0. For the NS 
box systems m — 1,2, ... , M, where the energy dependent M is the number of 'transverse modes'. For the NS disk 
systems \m\ < M, where M is the maximum of the angular momentum quantum number m for the Andreev states. 
We consider only the positive energy spectrum, therefore n > 0. Here we take the box geometry. The following 
expressions for disk geometry can be obtained straightforwardly by taking into account the fact that the set of m is 
different in this case. Thus M has a different meaning for box and disk geometries. We shall always indicate how 



to convert the box geometry results to the disk geometry case. To get Eq. (42), the phase term in Eq. ( |26| ) due to 
Andreev reflection, arccos(e/Ao) was approximated by tt/2, which is valid for e — > 0. Later a better approximation 
will be given by Taylor expanding the phase term, too. The density of states is then 
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M 



M 



P (e) = d ( e ~ £n - m ^ = / dsS [ £ 



n—O m—1 



n—0 m—1 



(43) 



We have seen in Sees. Ill A and MB that ?i$ m (0) is the classical action for the electron with Fermi energy moving 
in the N region between two subsequent bounces at the superconductor. Thus, s = hvF$' m (0) is the path length of 
the trajectories of the electron between two successive bounces at the NS interface. We now show that the above 
expression for p(s) can be rewritten in the same form as in (f4l|), which makes it possible to express the probability 
P(s) in terms of $ m (e). r-.r— . 

Applying the Poisson formulaEJ'Ea for the summation over to in (f43) we have 



n=0 



p(e) = E / dsSU- 



(n + flTTVF 



E 



M+1/2 



1/2 



dm 5 (s — hvF& m (0)) e 



i27rkm 



(44) 



Keeping only the non-oscillating term (fe = 0) in the sum over k, and then performing the integral over to we obtain 
the same Bohr-Sommerfeld approximation to the DOS as in Eq. (Mil). The probability P(s) is given by 



P(s) 



8 (M - m*) 9 (to* - 1) 



where the s-dependent to* satisfies 



Mhv F 



d<S>' m (0) \ 

dm. 



(45) 



(46) 



Q(x) is the Heaviside function. Note that the probability P(s) is normalized to 1, i.e. P(s)ds = 1. For disk 
geometry the numerator in (EjJ) has to be replaced by 0(M — |to*|)/2, and the necessary replacement in ( pl| ) is 
M — > 2M. Our main result is that the probability P(s) is expressed in terms of $ m (e), which has already been 
determined for box and disk geometries of the NS system (see Eqs. ( |30| ) and (A9)). 
Performing the integral over s in Eq. (f4l|) one finds 



PBs(e) = — s n (e)P{s n (e)), 



(47) 



where 



s n (e) 



Similarly, the counting function Nbs(s) — Jq de' pbs(s') can be easily found 



oc <>oo 
n=0 • /s «( £ ) 



(48) 



(49) 



Both in (^7|) and (^) we have to replace M by 2M in the case of NS disk systems as M has a different meaning 
in this case. The above results are valid only for small e. Numerical results show that a better approximation 
of the counting function can be obtained by Taylor expanding the phase shift related to the Andreev reflection: 
arccos(e/A ) ~ ir/2 — e/A . Then the discrete energy levels are given by 



Carrying out the same procedures as before one finds that the counting function N-bs(c) is still given by Eq. 
the replacement 



(50) 
after 



i s m ade in 



s n (e) -*■ s„(e) - Co (51) 
Here £0 = £( e = 0) = ?ivf/ 'Aq = 2Ep/ (kp^o) is the coherence length at the Fermi energy (see Eq. 
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In the following two subsections we shall calculate P(s) and the counting function -/Vbs(e) for NS box and disk 
systems. We shall see that the probability P(s), more precisely sP(s), is singular at some s s i ng for both NS systems. 
Then, if s S mg is known, one can determine from (|5|) the energies where P(s n (s)) — and consequently (see Eq. (pTij)) 
the density of states /Obs(e) — is singular: 

4^ = fe^A , (52) 

which is valid for all n for which Sn < Ao- Note that this expression can be applied to those normal systems 
attached to a superconductor for which sP(s) is singular. It is clear from this expression that these singularities are 
at equal distances. The analytical behavior of P(s), and thus, the existence of the singularities in the DOS for perfect 
NS interface inherently depends on the geometry of the isolated normal billiard. Once sP(s) is a singular function of 
s (because of the shape of the normal region), then the subsequent singularities in the DOS are at equal distances. 
We shall see that from the exact calculation of the energy eigenvalues for NS box and disk systems the positions 
of the singularities in the DOS agree very well with those given by the above derived approximate expression. We 
would like to emphasis that the second term in the denominator of Eq. ( |52"|) resulted from the Taylor expansion of 
arccos(e/A ) which is the phase shift due to the Andreev reflection. We shall see below that this extra term involving 
the coherence length £ in Eq. (|52| ) is necessary for getting good agreements between the exact and that of predicted 
by Eq. (p2[) for the position of the singularities. 



A. Singularities in DOS for Box 



In Sec. |III A| we have calculated $ m (e) for the NS box system with perfect interface. Inserting $ m (e) given by fl30| ) 
into Eqs. (}45|) and ( flo]) one finds 



poo 



Ad 2 



: Q(S - S„ 



(53) 



l-(f) 



where s m ; n 
modes, i.e. 



= 2d/y/l 

Mj-2> 1 we have s 



1/M 2 and M = k-pW/ir is the number of 'transverse modes'. For a large number of transverse 



2d. 



Notice that the same result can be found for P(s) by simple geometrical 
considerations^. P(s) is singular at s s ; ng = 2d, therefore the DOS is singular at the energies given by (|52|). The DOS 
can be calculated from Eq. ( [f7| ) in Bohr-Sommerfeld approximation. In Fig. ^the normalized DOS, pBs( e )/(2p ( ' N ' ) ) is 
plotted in Bohr-Sommerfeld approximation together with its slope for e — > 0. Here p' N ' = 2m/7i 2 ( J 4/47r) is the DOS 
of the normal box of area A — dW. 




FIG. 7. The Bohr-Sommerfeld approximation of the DOS for the box geometry — normalized by the DOS of a normal billiard 
of area A — dW (solid line) — and its slope for e — * (dashed line) as functions of e/Ao. The parameters are the same as in 
Fig. ti. With these parameters 2d/£o = 16.57. 



From (|48|) it is clear that the DOS for small e is dominated by the large s behavior of P(s). According to (|53j), 
P(s) — > 4d 2 s~ 3 in this case. Thus, 
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(54) 



where p (N) = (2m N /?i 2 ) A/ (4m) is the DOS for the isolated billiard of area A = dW and E T = M/(4irp iN) ) is the 
Thouless energy defined in Refs. [|fL5| This result is shown by the dashed line in Fig. [7| iJjbr small energies e the 
DOS is proportional to e in agreement with the findings of Melsen et al.BliJ and Ihra et alJl3. However, the slope of 
Pbs(s) /2p^ N ^ is less than what these authors found. The reason is that the geometries of the billiards they studied are 
slightly different from our box geometry. Namely, the width w of their superconducting lead is smaller than the side 
length of the billiard with which it is in contact. In these geometries, by quantum mechanical calculations, Ihra et 
al.ta found 3 pronounced peaks in the DOS approximately at s/Et — 3.1, 8.3, 13.4 (see Fig. 3. in that paper). They 
mention that similar peaks were observed for other parameter values of the billiards. If we assume that these peaks are 
indeed singularities in the DOS then fitting these data to Eq. ( |52j) one finds that P(s) is singular at s = s s i ng ~ I.ILt, 
where Lt = na 2 /w is the Thouless length used in their paper. The Thouless length is related to the mean escape time 
by T esc = Lt/vf for the billiard which is open along the superconducting lead. Since it is reasonable to expect P(s) 
to be high around s — Lt, our theory about the singularities in the DOS may explain the reason for the pronounced 
peaks in the DOS observed by Ihra and his coworkerstl They have calculated P(s) for another set of parameters of 
the billiard, and a peak can also be seen around s = Lt in Fig. 2. in their work. However, to confirm the existence 
of the singularity in P(s) one needs to calculate P(s) on a finer scale than in Fig. 2. of their paper. 

To avoid the errors of numerical differentiation, we shall compare the counting function obtained from the quantum 
mechanical calculations (see Eq. Jl9|)) with the counting function Nbs(z) i n Bohr-Sommerfeld approximation. Using 
( p9| ) and (|5^), and performing the integration, we obtain 

oo 

N BS (e) = Mj2f(sn(z)), (55) 

n=0 

where 

/(S) = { 1- v/l-4d 2 / S 2 , if s > Id, (56) 

and s n (s) is given in Eq. (|5l"|). The above derived counting function Nbs( s ) is plotted in Fig. || together with the 
exact counting function. 



son 




FIG. 8. The exact counting function N(e) (solid line) and Nbs(s) given in (^5|) (dashed line) as functions of e/Ao for the 
box geometry. The parameters are the same as in Fig. ti. 

One can see that the agreement is excellent at low energies (below the second cusp) but for higher energies the 
exact counting function is smaller than the approximated one. However, the positions of the singularities agree quite 
well in the two curves. Slight deviations are observed only for the last two cusps. The same was found for other 
parameter ranges. Thus, our expression (|52|) for the positions of the singularities in the DOS results in an excellent 
agreement with exact quantum mechanical calculations. Our result may highlight the origin of the singularities of the 
DOS. One can see that P(s) is singular at s = 2d, which corresponds to the classical orbits when the electron enters 
the box parallel to the x axis and is reflected back at the wall. The enhanced return probabilities for these orbits are 
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quite obvious. These are the orbits which result in the singularities of the DOS. It is also clear from this expression 
that the singularities are at equal distances. 

The expression for the positions of the singularities given by Eq. ( |52| ) can be applied to other NS systems, e.g. 
that studied by de Gennes and Saint- JamesQ. They considered an NS system in which a normal film of width a is in 
contact with a semi-infinite superconductor. This geometry of the NS system is similar to our NS box system shown 
in Fig. 0a. The authors also found singularities in the DOS for parameter values 2a/ = 2.5,8,20. Using Fig. 1 of 
their paper we measured the positions of the singularities in the DOS, labeled them as n = 0, 1, 2, • • •, and plotted 
them in Fig. ^ (circles and squares) against n. For 2et/£o = 2.5 they found only one singularity which is not shown in 
Fig. |^. We can expect that the probability P(s) for this system, similarly to our box geometry, is singular at s = 2a. 
Thus, assuming that s s i ng = 2a we calculated the positions of the singularities from Eq. ( |52| ) and also plotted them 
in Fig. |9|. For clarity, these points are connected by lines. 



.2^=20 
- BS 
■ 2a/^=8 
BS 



I 1 1 1 1 1 1 1 1 

-1 1 2 3 4 5 6 

n 

FIG. 9. The positions (in units of Ao) of the singularities of the DOS obtained from Fig. 1 in the work of de Gennes and 
Saint- JameaD (circles and squares) and from our expression ( |52| ) (for clarity, the points are connected by lines). 

It can be seen from the figure that the agreement is excellent. Even for 2a/£o = 2.5 (when there is only one 
singularity) the position of the singularity agrees very well with that found from our expression (|52]). We note that 
to achieve such a good agreement the Taylor expansion of the phase shift arccos(e/Ao) in terms of e was necessary. 
Therefore, to understand the nature of the singularities one needs to take into account not only the singular behavior 
of P(s) which depends on the geometry of the billiard but a more physically related quantity, namely the coherence 
length of the superconductor. .— . 

Finally, we mention that Lodder and NazarovL3 also found singularities in the DOS for different shapes of Andreev 
billiards. Our results suggest that these singularities are related to some special classical orbits which are characteristic 
of the specific geometry. Sipr et al.EJ studied SNS systems taking fully into account the motion parallel to the infinite 
interface. Their results also show singularities in the DOS. We believe that these singularities are also related to 
some special classical orbits which give rise to a singular behavior of the probability P(s), similarly to the NS systems 
studied in this paper. 

We also calculated the energy levels in the case of non-perfect NS interface. Solving the secular equation ( |l9| ) 
numerically, the obtained counting function N(e) is shown in Fig. |l0[ We have used the same parameters for describing 

the mismatch in the Fermi energies and the effective masses as in Ref. |3^, namely ru = /kp\ r v = /v^\ The 
strength of the tunnel barrier is given by Z. 
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FIG. 10. The counting function N(e) as a function of e/Ao (solid line) for the box geometry when the interface is not perfect. 



The parameters are: d/W 



3, k^W/n 



87.9, A /E 



(N) 



0.02, Tk = 0.007, r v =0.1 and with no tunnel barrier, i.e. Z = 0. 



The Weyl formula in leading order (see text) is given by the dashed line. 



It can be seen from the figure that for a non-perfect NS interface the cusps in the counting function are 'washed out'. 
It is related to the fact that the normal reflections at the interface are enhanced. The same is true for Z ^ 0. In Sec. 
Ill A we have found that in leading order the Weyl formula is given by N(e) ~ 2e p( N \ where p( N ) = 2m,/h 2 (A/ An) is 
the density of states for the isolated N region of area A = Wd. This is shown by the dashed line in Fig. [l(]. One can 
see from the figure that the exact counting function can be approximated by the above Weyl formula. This implies 
that for non-perfect NS interfaces the effect due to the Andreev reflection is less significant, and the system behaves 
as a normal metal regarding the energy levels. 



B. Singularities in DOS for disk 



The maximum angular momentum quantum number, A/as for Andreev states has been determined in Sec. [ill B| 
after Eq. (34). To have a better physical picture of these states we now give an estimation for Mas based on 
semiclassical considerations. It is easy to see that the angular momentum (in units of ft) of the rays shown in Fig. ^a 
is m — fcpi?N shi(ev/2), where a is the angle between the two rays at the outer circle. Rays that are tangent to the 
inner circle have maximal angular momentum, in which case sin(a/2) = Rs/R^. These are the trajectories which 
still reach the superconductor so that they belong to the Andreev states. Thus, the maximal angular momentum m 
for the Andreev states is k^Rg. This is a good estimate for Mas given after Eq. (|34|). 

In subsection IIIB we have seen that for Rg ^> £(e) and for perfect NS interface the energy levels for NS disk 



systems can be classified into Andreev and whispering gallery states to a good approximation. We first consider the 
Andreev states. Inserting <& m (e) given by Eq. (|A9[) into Eqs. p5|) and (46), one finds after some algebra 



P(s) 



1 



4 s 2 R s x /p 7s 

° V 1 max/ J 



s 2 . 

mm 



[s 2 



:0( Sn 



s)6(s 



0, 



(57) 



where s m m = 2(i?N — Rs) an d s max = 2y/ R 2 ^ — i? 2 ,. The two Heaviside functions come from the factor 8(M— |m*|)/2 
(see the text after Eq. ([46|)) when we change to the variable s. Here s m i n and s max correspond semiclassically to the 
smallest and the largest path lengths of the orbits related to the Andreev states. Obviously, s min equals to the path 
length of those orbits for which the trajectory of the electrons between two successive bounces at the NS interface is 
along the radius of the two concentric circles. The trajectories of the electrons corresponding to the maximal path 
length for Andreev states are those which are tangent to the inner circle of radius Rg. The probability P(s) is zero 
for s m in > s > Smax and normalized to one. Note that in deriving (57), the solution of Eq. ( (46| ) is twofold for m* 
resulting in aia-extra factor 2 in P(s). Notice that the same result can be found for P(s) by simple geometrical 
considerations^ 2 ! One can see that P(s) is singular at s = s m i n . The contributions of the Andreev states to the DOS 
in Bohr-Sommerfeld approximation is given by Eq. (|47|), which reads 



2ki? Rs 

Pad (s) = s n (e)P(s n (e)), 



(58) 



n=0 
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where P(s) is given by (57) and s„(e) by (pl[) . Since there is a maximum path length s for Andreev states, the sum 
over n in ( [47] ) is indeed finite. 

Next, we consider the contribution of the whispering gallery states to the DOS. After differentiations with respect 
to e in Eq. ( |37| ) and replacing the summations over m by integrals the calculations can be carried out analytically 
(here we used M wgs sa k e Rg and neglected the higher order terms in e/Ep). It is found that the DOS is constant in e 
and is given by 



(k F R N y 
2E F 




arcsm ■ 



Rs 



(59) 



Thus, the total DOS for the NS disk system in Bohr-Sommerfeld approximation is p wgs + pad(s)- Since P(s) is a 
singular function, the DOS becomes singular at the energies given by ( |5^ ) with s s i ng = 2(i?N — -Rs)- In Fig- [ll] we 
plotted the normalized DOS, pbs( £ )/(2/5 ( - N ' ) ) in Bohr-Sommerfeld approximation, where = 2m/?i 2 (A/47r) is the 
DOS of the normal circular annular billiard of area A = (R^ — R^)tt. 



FIG. 11. The normalized DOS (see text) in Bohr-Sommerfeld approximation as a function of e/Ao for the disk geometry. 
The parameters are the same as in Fig. H. With these parameters s m i n /£o = 15.0. 



One can see that the DOS does not tend to zero as e — * 0, contrary to the box geometry case. The non-vanishing 
DOS is due to the constant contribution of the whispering gallery states given by Eq. (p9h. 



Again, one can determine the counting function A^es(e) by integrating the DOS. Using Eqs. (49) and (|57| ) the 
integral can be performed analytically, and 



is obtained, where 



Nbs(s) = Pwgse + 2k F R s J2 /Me)) ©(« 



s„(e)) 



n=0 



(60) 



m = 



2 ][ 



s*\\s A —s* 



if S ^ Smini 

if Smin < 8 < S r 



(61) 



and s n (e) is given by (pl| ). Note that the infinite sum over n in ( J60| ) is indeed a finite one due to the Heaviside 
function. Using Eqs. (|19[)-(p0[) we calculated the exact energy levels for the NS disk system. In Fig. [l^ the exact 
counting function N(e) and Abs(^) are plotted for perfect interface (no mismatch and zero tunnel barrier). 
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0.4 0.6 

FIG. 12. The exact counting function N(s) (solid line) and Nbs(s) given in ( |6o| (dashed line) as functions of e/Ao for the 
disk geometry. The parameters are the same as in Fig. o. Then s m i n /£o = 12.5. 



One can see that the agreement between the exact counting function and that obtained from the Bohr-Sommerfeld 
approximation is excellent for all energies below the gap. To see the agreement, the details of Fig. |l2|- up to the first 
two cusps - are enlarged in Fig. |l^. 




FIG. 13. Enlarged portion of Fig. 

Similarly, very good agreements were found for other parameter ranges of the NS disk systems with perfect interface. 

In Fig. I one can see that the energy levels (for each radial quantum number n) as functions of m are approximately 
constant as m — ► 0. Thus, the DOS should diverge because it is proportional to the reciprocal of the derivative of e mn 
with respect to m (assuming that m is a continuous variable). The positions of the singularities of the DOS coincide 
with the energies e mn at m = for all possible n. The angular momentum quantum number m — semiclassically 
corresponds to the radial orbits of the electron between the inner and outer circles. Applying the Bohr-Sommerfeld 
quantization rule for the electrons and the Andreev reflected holes along the radius one finds 

[k e (e) — kh(e)} (i?N — Rs) = n7T + arccos(e/A ). (62) 

Note that this equation can also be derived from the quantization condition given in Eq. (26). The solutions of Eq. 
( |u2| ) agree very well with the energies in Fig. [| at m — 0. Therefore, the singularities of the DOS are related to the 
classical orbits along the radial directions. 

We also calculated the energy levels in the case of non-perfect NS interface. Solving the secular equation jl9| ) 
numerically, the counting function N{e) is shown in Fig. 14. 
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Exact 

Leading Weyl 




FIG. 14. The exact counting function N(e) (solid line) and the leading term of the Weyl formula N(e) (dashed line) as 
functions of e/Ao for the NS disk geometry when the interface is not perfect. The parameters are Rs/Rn = 7/10, Rs = 350, 
A /£4 N) = 0.1, r k = 0.01, r v = 0.1, and tunnel barrier, Z = 0. 

It can be seen from the figure that for non-perfect NS interface the cusps in the counting function are 'washed out' 
similarly to the case of NS box systems. This is related to the fact that the normal reflections at the interface are 
enh anced. Our numerical results show that the cusps also disappear for Z ^ 0. Using the Weyl formula derived in 



Sec. Ill B we found that N(e) = 2p( N ) e in leading order, where p' N ' = (2m/?i 2 )(i?^ — i?g) /4 is the density of states 
for normal annular billiard with inner radius Rs and outer radius R^. In Fig. n4]this is drawn by a dashed line. One 
can see that this Weyl formula is a good approximation for the exact counting function in the case of non-perfect 
interfaces. The result shows that the role of normal reflections dominates over Andreev reflections, and thus the Weyl 
formula agrees approximately with that for the circular annular billiard without superconducting region. 



V. 'PHASE DIAGRAM' FOR NS DISK SYSTEMS 



In our studies of the energy levels of the NS disk systems we have assumed so far that the energy dependent 
coherence length £(e) = hvp/ yAg — e 2 is much smaller than Rs (see the text after Eq. flA8|)). In this section we 
shall discuss the case when ^> Rs, i.e. when e — > Ao. The imaginary part of q e Rs ~ kpRs + iRs/£,( e ) is then 
negligible, and k e Rs ~ q e Rs- Thus, in the determinant (e) given in Eq. ( |23| ) the arguments of the Bessel functions 
are almost equal. We now consider the perfect NS interface (Z — and ton = to§). Subtracting the second column 
from the first one in the determinant, the ratio J m (k e R^) /Y m {k e R^) can be factored out. It turns out th at the 
remaining determinant is non-zero for e < Ao- Thus, similarly to the whispering gallery states discussed in Sec. Ill B , 
it holds to a very good approximation that Dm \s) (as a function of e) has zeros whenever J m (k e R^) = for all to. 
This corresponds to the energy levels of a circular billiard of radius i?N for electron-like quasiparticles. The same is 
true for the hole-like quasiparticles, i.e. their energy levels are the solutions of JmikhRn) = for all to. In summary, 
we found that for £(e)/i?s S> 1 the NS disk system behaves as a normal full disk, while in the opposite case the energy 
levels of the sy stem have a 'mixed phase' character, in which Andreev states and whispering gallery states coexist 
(see Sec. IIIB). We shall call the first case a 'full disk phase' (FD). According to our numerical results, it is safe to 
say that the MP and FD can be distinguished very well by the ratio £,(e)/Rs- We found mixed phase for £(e)/i?s < 2 
and FD phase for £(e)/R s > 2. 

As it is known, the Andreev approximation is valid only for Aq/Ep < 1. If we assume some critical value, say 
Ao/E-p =0.1, then for fcp£o = 2£ , p/Ao < 20 the Andreev approximation is not valid and normal reflection at the 
NS interface dominates over Andreev reflection. The system looks like an annular billiard (without superconductor 
region). We shall call this case an 'annular phase' (AP). As we have seen in the Fig. [l4|, the non-perfect interface also 
results in an enhanced normal reflection and the character of the energy levels again corresponds to a circular annular 
billiard (annular phase). 

We calculated the exact energy levels for different parameters and compared with those obtained for the full disk 
and for circular annular disk. In this way, we can classify the NS disk systems into the MP, FD and AP phase regions. 
The crossovers between the three 'phases' are not sharp and depend on the energy, too. It turns out that the three 
phases can be characterized by two parameters, kpRs-, kp^(s)- These are the relevant parameters for the NS disk 
systems and each pair corresponds to one of the phases. Thus, a so-called phase diagram can be given for the different 
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phases. The boundaries of these phases are not so sharp, however. Far from the boundaries the given phase becomes 
dominating. In Fig. [l^ the three different phases are shown. 




^(8) = 2R S 



MP 



AP 




FIG. 15. Schematic phase diagram for the NS disk systems. The mixed phase (MP) is bounded by the lines £;f£o ~ 20 and 
£(e) w 2R S . The full disk phase (FD) is bounded by the line k F R s = 0, the line fc F £o ~ 20 and the line £(e) « 2R S . The 
annular phase (AP) is below the line £;f£o ~ 20. 



We investigated the energy spectrum of NS box and disk systems using the BdG equation. Matching the wave 
functions at the NS interface we derived a secular equation whose solutions give the energy spectrum of the system. 
The mismatch in the Fermi wave numbers and the effective masses of the normal system and the superconductor, 
as well as the tunnel barrier is included in the calculation. Rewriting the secular equation a simple quantization 
condition was derived. Equation (^6|) is a central result and serves as a starting point for further theoretical studies 
in this paper. Using this quantization condition (a) we derived for both NS systems with perfect interfaces a Weyl 
formula for the counting function, (b) by re-deriving the commonly used Bohr-Sommerfeld approximation to the DOS 
we presented an explicit expression for the probability P(s) in terms of the classical action of the electron moving in 
the N region between two successive bounces at the NS interface. The numerically exact calculations are in very good 
agreement with our Weyl formula and the Bohr-Sommerfeld approximation. The singularities in the DOS obtained 
from our exact calculations can be successfully explained by using the correct probability function P(s). P(s) is a 
singular function of s both for NS box and disk systems. According to our theory concerning the singularities, this 
implies that the singularities of the DOS are located at equal distances from each other. A simple formula is given far 
the location of the singularities. We demonstrated that it can be applied for the system studied by de Gennes et al.Lj, 
and it may as well give the reason for the significant peaks observed by Ihra et al.lij in their numerical calculations. 
We show that the singularities in P(s) are related to some special classical trajectories of the electron moving in the 
N region hitting the NS interface. 

In the case of NS disk systems the DOS is constant as the energy gets close to the Fermi level. This is, at first 
sight, in contrast to the common belief that for integrable billiards the DOS should go to zero as the energy tends the 
Fermi level. We pointed out that in this NS disk system the whispering gallery states give the non-vanishing DOS at 
the Fermi level. However, the Andreev states still have no contribution to the DOS. In the disk system the Andreev 
and the whispering gallery states coexist (so-called mixed phase). Depending on the geometry of the disk and the 
material parameters, we sketch a kind of phase diagram to classify the energy spectra into three classes: (i) mixed 
phase, (ii) full disk, (iii) annular disk. If the coherence length is much larger than the diameter of the superconducting 
region, the spectrum corresponds to the full disk phase. The electrons travel thorough the S region without Andreev 
reflection. When the Andreev approximation fails, the normal reflection is enhanced at the NS interface and the 
system behaves as an annular circular billiard. This is also the case when there is a tunnel barrier at the NS interface 
or a mismatch between the effective masses and the Fermi wave numbers of the normal and superconducting regions. 

An interesting extension of the study presented in this paper may be the investigation of the role that the geometry 
of the normal billiard plays in the singularities of the DOS. Under what conditions do such singularities exist? Does 
the energy spectrum at higher energies show some special structure in chaotic billiards? How can the Weyl formula 
(derived for NS box and disk systems) be extended for other integrable and chaotic billiards? Our quantization 
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method can also be applied for NS disk in the presence of a magnetic field. In this case, besides Andreev states and 
whispering gallery states (more specifically edge states) one has to take into account the Landau levels that appear 
in the spectrum in the presence of strong magnetic fields. These give rise to a much more complex energy spectrum 
than a zero magnetic field. Another relevant question is the study of the counting function in these systems. Further 
work along these lines is in progress. 
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APPENDIX A: CALCULATION OF THE PHASE $ m (e) 

It is clear for symmetry reasons that the eigenstates with quantum numbers m and — m are degenerate (except for 

m = 0), thus it is enough to consider the states m > 0. In the Andreev approximation (Ao/E-p <C 1) we take k$ « q$ 

in places where they are multiplied by the Bessel functions in Dm (e) but in the arguments of the Bessel functions 
we keep them to be different. As we shall see, different methods are necessary to approximate the determinant 

D { rn>(e) for different ranges of m. The Debye asymptotic expansion!^ will be used to approximate the determinant 
for |m| < k e Rs — \/k e Rg. In this expansion the Bessel functions with the real argument x for < m < x — <jx~ are 
approximated by 



1 



7r y/x sin Q 



2 1 



7T 



Jm( x ) ~ V cos x (sinQ — QcosQ) — — , (Al) 



7r y/x sin Q 



TT 



Y m (x) w \ . sin x(sinQ - QcosQ) - - , (A2) 



4 



where cosQ = m/x. For \m\ > k e Rs + \fk e R§ an analogous expansjonii can be used. In the intermediate range 
|m — k e Rs\ < \/k e Rg one may apply uniform asymptotic expsjisionsE-J for the Bessel functions. For Bessel functions 
of a complex argument z with fixed m the approximate forma 

J m (z) ~ \ — cos ( z mir ] , for z — > oo (A3) 

V nz V 2 4 / 



will be used. We first calculate the eigenphase $ m (e) for \m\ < k e Rs — \/k e Rs. Using Eqs. ( Al )-( |A3| ) for the Bessel 
functions in Eqs. (|3|) and (p5|), one finds 

1 i p -2i[f„ l (fc c fl N )-^ m (fc c fl s )+/3*] 

p i*m{e) _ 1 + e 2 l [tf m (fee_R N )- 1 9 m (fe _R s )] p i(f3+0») ( aa\ 

I + e 2i[ 1 ? m (fc e fl N )-fl m (fc« ! _Rs)+/3] C ' 

where 

$m( x ) — V x 2 — m 2 — |m| arccos (A5) 

1 7T 

(3 = q e R s - -mw - -. (A6) 

The absolute values of m in the above expressions are the consequence of degeneracy of the eigenstates m and — m. 
Note that j3 is a complex number. Thus, for Ao <C E-p and below the gap 77 <§C 1 , Eq. ( |l2|) yields 

1 tt Rs 

(3 = k F R s - -rrnr - - (A7) 
where the coherence length ^(e) is definecfi as 
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i (e) = 7=^== = rTT A8 

Here up is the Fermi velocity. We now assume that i?s The opposite case will be discussed in Sec. [v|. Then, 

the first factor in Eq. (A4) is approximately equal to 1, and we obtain 



$ m (e) = 2 [# m (k e R N ) - ^m(fce-Rs)] + 2fc F i? s -rrm-1. (A9) 



We now turn to the case \m\ > k e Rg + \/k e Rs- In this case, according to the analogous Debye asymptotic expansion 
of the Bessel functions, one finds that J m {k e Rs)ca,n be neglected (it is exponentially small) compared to the second 
term of the 1,1 element of the determinant in (p3|). Similarly, J' m (k e Rs) is also negligible in the determinant. Thus, 
the ratio J m {k e R^{) /Y m {k e R^) can be taken out from the determinant. It turns out that the remaining determinant 

is not equal to zero for s < Aq. Thus, to a very good approximation, J)t' (e) (as a function of e) has zeros whenever 
J m {k e R^) = 0. This corresponds to the energy levels of a circular billiard of radius R^ for electron- like quasiparticles. 

Note that now $ m (e) cannot be calculated from Eq. (p5|), since Dm(s) ~ for \m\ > k e Rs + v / fc e i?s- In this case we 
shall use a different method to determine the contributions to the smooth part of the counting function. 

To approximate D$(s) in the intermediate range |m — k e Rs\ < \/k e Rs, one may apply uniform asymptotic 
expansions for Bessel functions. Semielassically, these states are related to the diffraction of the electron in the 
penumbra of a circular annulus billiaroS. As a first estimate of the smooth part of the counting function, the energy 
levels in the intermediate range will be taken into account by the corresponding energy levels of the circular billiard of 
radius i?N using Debye's asymptotic expansion. As we shall see, this is a good approximation since the intermediate 
range of m is rather narrow compared to the two other ranges discussed above. Using the uniform expansions of 
the Bessel functions in our derivation of the Weyl formula for the NS disk system, it is straightforward (although 
algebraically rather tedious) to include the intermediate range more accurately. 

Similarly, the method of approximations outlined above can be used to approximate (e) in the secular equation 



IS). One finds that <E> m (— e) is still given by (A9) after replacing k e by the wave number kh(e) = k e (—e) of 



the hole-like quasiparticles in the N region. However, in this case the approximation for Q m (— e) is valid only for 



|m| < khRs — ykhRs- Finally, from Eq. flA9| ) one obtains 

$m(e) - $m(-e) = 2 [$ m (k e R N ) ~ # m (k e Rs)] - 2 [0 m (k h Rs) - $ m (k h R s )} (A10) 



for |m| < khRs — vfaiRs- Note that the last three terms in (JA9J) cancel out in the difference $ m (e) — $ m (— e). 

In a similar way, it can be seen that the zeros of Dm \e) in Eq. ( |l9| ) coincide with the zeros of J m (khRN) for 
\m\ > khRs + v 7 khRs to a good approximation. Thus, the secular equation (|l^), determining the energy levels of the 
NS disk system, reduces to 

J m (k e R N )J m (k h R N )=0 (All) 



for |m| > k e Rs + \/k e Rs. The energy levels are the same as those of the disk of radius i?N for electron/hole- like 
quasiparticles. Semielassically, the states with angular momentum quantum number |m| > fc e i?s + ykeRs are called 
whispering gallery modes. The wave functions for these states at r = Rs are negligible, and so no Andreev reflections 
take place. This is why the superconducting pair potential A does not appear in ( All] ). 
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